function f = custom_fpdf(x,v,w)
% custom_fpdf - PDF of F (Fisher-Snedecor) distribution
%
% FORMAT:       f = custom_fpdf(x, v, w)
%
% Input fields:
%
%       x           F-variate
%       v           parameter 1 / numerator d.f. (v > 0)
%       w           parameter 2 / denominator d.f. (w > 0)
%
% Output fields:
%
%       f           PDF of F distribution with v,w d.f. at points x
%
% Note: this function has been copied from the SPM2 package published
%       by the Wellcome Department, go here for details:
%       http://www.fil.ion.ucl.ac.uk/spm/

% Definition:
%-----------------------------------------------------------------------
% The PDF of the F-distribution with degrees of freedom v & w, defined
% for positive integer degrees of freedom v>0 & w>0, and for x in
% [0,Inf) by: (See Evans et al., Ch16)
%
%             gamma((v+w)/2)  * (v/w)^(v/2) x^(v/2-1)
%    f(x) = --------------------------------------------
%           gamma(v/2)*gamma(w/2) * (1+(v/w)x)^((v+w)/2)
%
% Variate relationships: (Evans et al., Ch16 & 37)
%-----------------------------------------------------------------------
% The square of a Student's t variate with w degrees of freedom is
% distributed as an F-distribution with [1,w] degrees of freedom.
%
% For X an F-variate with v,w degrees of freedom, w/(w+v*X^2) has
% distributed related to a Beta random variable with shape parameters
% w/2 & v/2.
%
% Algorithm:
%-----------------------------------------------------------------------
% Direct computation using the beta function for
%       gamma(v/2)*gamma(w/2) / gamma((v+w)/2)  =  beta(v/2,w/2)
%
% References:
%-----------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
   ~isnumeric(x) || ...
   ~isnumeric(v) || ...
   ~isnumeric(w) || ...
    isempty(v) || ...
    isempty(w) || ...
    any(isnan(x(:)) | x(:) < 0) || ...
    any(isinf(v(:)) | isnan(v(:)) | v(:) <= 0) || ...
    any(isinf(w(:)) | isnan(w(:)) | w(:) <= 0)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(x)
    f = [];
    return;
end
osize = size(x);
if numel(x) == 1
    if numel(v) > 1
        osize = size(v);
    elseif numel(w) > 1
        osize = size(w);
    end
end
onum = prod(osize);
if ~isa(x, 'double')
    x = double(x(:));
else
    x = x(:);
end
if ~isa(v, 'double')
    v = double(v(:));
else
    v = v(:);
end
if ~isa(w, 'double')
    w = double(w(:));
else
    w = w(:);
end
if onum > 1
    if numel(x) == 1
        x = repmat(x, [onum, 1]);
    end
    if numel(v) == 1
        v = repmat(v, [onum, 1]);
    end
    if numel(w) == 1
        w = repmat(w, [onum, 1]);
    end
end
if numel(x) ~= onum || ...
    numel(v) ~= onum || ...
    numel(w) ~= onum
    error( ...
        'BVQXtools:SizeMismatch', ...
        'Invalid number of elements in x, v, or w.' ...
    );
end

% compute and reshape
f = reshape((v ./ w) .^ (v ./ 2) .* x .^ (v ./ 2 - 1) ./ ...
	(1 + (v ./ w) .* x) .^ ((v + w) ./ 2) ./ beta(v ./ 2, w ./ 2), osize);
